Covid-19 Heatmap & Population_Density

利用近期新冠疫情於台灣之確診者足跡製作熱點圖(核密度法,Kernel Density Estimation），並與台北市人口密度進行套疊，以觀察確診者足跡熱點空間分布與台北市人口密度的關係。

import pandas as pd 
import geopandas as gpd
from pyproj import CRS
import requests
import geojson
import folium

村里人口數

STEP1 使用臺北市政府民政局統計資料 臺北市每月各里人口數及戶數

STEP2 新建一個TownVillage欄位，為區域和村里的文字合併(作為後續關聯使用)

# STEP1 加入統計資料
taipei_pop = pd.read_csv(r"E:\Py_NTU 341\Homework\Data Folder\110年4月台北市各里人口數戶數.csv",encoding='big5')
taipei_pop.head()

	區域	村里	村里數_現有門牌	村里數_戶籍登記	鄰數_現有門牌	鄰數_戶籍登記	戶數	人口數-合計	人口數-男	人口數-女
0	松山區	莊敬里	1	1	27	27	2056	5279	2548	2731
1	松山區	東榮里	1	1	33	32	3071	8076	3773	4303
2	松山區	三民里	1	1	32	32	2875	6760	3086	3674
3	松山區	新益里	1	1	18	18	1843	4569	2174	2395
4	松山區	富錦里	1	1	23	23	2125	5181	2369	2812

# STEP2 新建一個TownVillage欄位，為區域和村里的文字合併(作為後續關聯使用)

taipei_pop["t.v"] = taipei_pop["區域"] + taipei_pop["村里"]
#taipei_pop.head()

pop = taipei_pop[["t.v","人口數-合計"] ]

pop = pop.rename(columns={'人口數-合計': 'pop202104'})
pop.head()

	t.v	pop202104
0	松山區莊敬里	5279
1	松山區東榮里	8076
2	松山區三民里	6760
3	松山區新益里	4569
4	松山區富錦里	5181

台北市村里圖資整理

STEP1 找到全台村里界圖 村里界圖(TWD97_121分帶)

STEP2 選出台北市之範圍

STEP3 計算村里面積(計算人口密度用)

STEP4 新建一個TownVillage欄位，為區域和村里的文字合併(作為後續關聯使用)

STEP5 計算人口密度

# STEP1  選出台北市之範圍
tw_village = r"E:\Py_NTU 341\Homework\Data Folder\VILLAGE_MOI_121_202104\VILLAGE_MOI_121_1100415.shp"
tw_village = gpd.read_file(tw_village,encoding = "UTF-8")
#tw_village.head()

# STEP2 選出台北市之範圍 
tp_village = tw_village.loc[tw_village['COUNTYNAME']=='臺北市']
tp_village.head()

	VILLCODE	COUNTYNAME	TOWNNAME	VILLNAME	VILLENG	COUNTYID	COUNTYCODE	TOWNID	TOWNCODE	NOTE	geometry
6952	63000080031	臺北市	文山區	樟新里	Zhangxin Vil.	A	63000	A11	63000080	None	POLYGON ((306484.609 2763701.609, 306442.073 2...
6953	63000080037	臺北市	文山區	老泉里	Laoquan Vil.	A	63000	A11	63000080	None	POLYGON ((307698.415 2763909.981, 307784.205 2...
6954	63000080032	臺北市	文山區	樟腳里	Zhangjiao Vil.	A	63000	A11	63000080	None	POLYGON ((307451.461 2763912.590, 307424.695 2...
6955	63000080041	臺北市	文山區	樟文里	Zhangwen Vil.	A	63000	A11	63000080	None	POLYGON ((306530.118 2763845.236, 306483.024 2...
6956	63000080043	臺北市	文山區	樟樹里	Zhangshu Vil.	A	63000	A11	63000080	None	POLYGON ((306790.748 2763911.814, 306732.793 2...

# STEP3 計算村里面積(計算人口密度用)

tp_village2 = tp_village.copy()   #避免py:853: SettingWithCopyWarning
tp_village2['Area'] = tp_village.area/10000  
tp_village2.head(2)

	VILLCODE	COUNTYNAME	TOWNNAME	VILLNAME	VILLENG	COUNTYID	COUNTYCODE	TOWNID	TOWNCODE	NOTE	geometry	Area
6952	63000080031	臺北市	文山區	樟新里	Zhangxin Vil.	A	63000	A11	63000080	None	POLYGON ((306484.609 2763701.609, 306442.073 2...	19.511016
6953	63000080037	臺北市	文山區	老泉里	Laoquan Vil.	A	63000	A11	63000080	None	POLYGON ((307698.415 2763909.981, 307784.205 2...	377.403620

#做 Town + village

tp_village3 = tp_village2.copy()
tp_village3["t.v"] = tp_village["TOWNNAME"] + tp_village["VILLNAME"]
tp_village3.head(2)
#village = tp_village[["COUNTYNAME","TOWNNAME","t.v","geometry","Area"]]
#village.head()

	VILLCODE	COUNTYNAME	TOWNNAME	VILLNAME	VILLENG	COUNTYID	COUNTYCODE	TOWNID	TOWNCODE	NOTE	geometry	Area	t.v
6952	63000080031	臺北市	文山區	樟新里	Zhangxin Vil.	A	63000	A11	63000080	None	POLYGON ((306484.609 2763701.609, 306442.073 2...	19.511016	文山區樟新里
6953	63000080037	臺北市	文山區	老泉里	Laoquan Vil.	A	63000	A11	63000080	None	POLYGON ((307698.415 2763909.981, 307784.205 2...	377.403620	文山區老泉里

# STEP4 新建一個TownVillage欄位，為區域和村里的文字合併(作為後續關聯使用)

tp_density = tp_village3.merge(pop,on="t.v")    #.fillna(0) 補0
tp_density.head(2)

	VILLCODE	COUNTYNAME	TOWNNAME	VILLNAME	VILLENG	COUNTYID	COUNTYCODE	TOWNID	TOWNCODE	NOTE	geometry	Area	t.v	pop202104
0	63000080031	臺北市	文山區	樟新里	Zhangxin Vil.	A	63000	A11	63000080	None	POLYGON ((306484.609 2763701.609, 306442.073 2...	19.511016	文山區樟新里	5662
1	63000080037	臺北市	文山區	老泉里	Laoquan Vil.	A	63000	A11	63000080	None	POLYGON ((307698.415 2763909.981, 307784.205 2...	377.403620	文山區老泉里	912

# STEP5 計算人口密度
tp_density["density"] = tp_join["pop202104"] / tp_join["Area"]
tp_density.head(2)

	VILLCODE	COUNTYNAME	TOWNNAME	VILLNAME	VILLENG	COUNTYID	COUNTYCODE	TOWNID	TOWNCODE	NOTE	geometry	Area	t.v	pop202104	density
0	63000080031	臺北市	文山區	樟新里	Zhangxin Vil.	A	63000	A11	63000080	None	POLYGON ((306484.609 2763701.609, 306442.073 2...	19.511016	文山區樟新里	5662	290.195032
1	63000080037	臺北市	文山區	老泉里	Laoquan Vil.	A	63000	A11	63000080	None	POLYGON ((307698.415 2763909.981, 307784.205 2...	377.403620	文山區老泉里	912	2.416511

#查看資料分布
plt = tp_density["density"].plot(figsize=(10,5))
plt.grid("on", axis="y")
plt.grid("on", axis="x")
plt.set_title("density")
plt.set_xlabel("VILLNAME")

Text(0.5, 0, 'VILLNAME')

#敘述性統計
tp_density["density"].describe()

count     455.000000
mean      284.789686
std       188.315154
min         0.545540
25%       141.283916
50%       280.230794
75%       403.193070
max      1112.840866
Name: density, dtype: float64

#展圖
tp_density.plot(figsize=(6,6),column='pop202104', scheme="natural_breaks",cmap='Greens', k=10)  #legend=True
tp_density.plot(figsize=(6,6),column='density', scheme="natural_breaks",cmap='Greens', k=10)  #legend=True

<AxesSubplot:>

製作covid案例熱點圖

covid19 = r"E:\Py_NTU 341\Homework\Data Folder\0516covid19_point.shp"
covid19 = gpd.read_file(covid19,encoding = "big-5")
#covid19.head(3)

covid19[["Name","geometry"]]

	Name	geometry
0	4/28 桃園市某國民小學	POINT Z (121.28966 25.04654 0.00000)
1	5/1Pilot in Cafe	POINT Z (121.20416 25.00696 0.00000)
2	5/1胡同彭家老舖 新疆拉麵	POINT Z (121.29158 25.04643 0.00000)
3	5/1華泰名品城	POINT Z (121.21327 25.01412 0.00000)
4	5/3 台北福華大飯店	POINT Z (121.54297 25.03747 0.00000)
...	...	...
272	教育部永續校園_新北市秀山國民小學	POINT Z (121.52114 24.99516 0.00000)
273	5/7 患者自行前往 診所就醫	POINT Z (121.52060 25.00017 0.00000)
274	三水街	POINT Z (121.50136 25.03608 0.00000)
275	西園路一段	POINT Z (121.49945 25.03709 0.00000)
276	雙和醫院－患者自行前往	POINT Z (121.49381 24.99130 0.00000)

277 rows × 2 columns

# Get x and y coordinates for each point
covid19["x"] = covid19["geometry"].apply(lambda geom: geom.x)
covid19["y"] = covid19["geometry"].apply(lambda geom: geom.y)

# Create a list of coordinate pairs
locations = list(zip(covid19["y"], covid19["x"]))

# 創建Choropleth map階級區分圖所需的Geo-id（每行必須具有唯一的標識符）
tp_density['geoid'] = tp_density.index.astype(str)

from folium.plugins import HeatMap


m = folium.Map(location=[25.04, 121.53], tiles = 'cartodbpositron', zoom_start=12, control_scale=True)


folium.Choropleth(geo_data = tp_density, 
                  data = tp_density,
                  name='Density',
                  columns=['geoid','density'],              #資料欄位
                  key_on='feature.id', 
                  fill_color='Greens',                      #Choropleth map色塊
                  line_color='white',                       #Choropleth map線色
                  line_weight=0,                            #Choropleth map線粗                   
                  legend_name= 'Density').add_to(m)         #圖例名稱
                        
    
#工具提示->顯示每個功能的值
folium.features.GeoJson(tp_density,
                        name='Density label',
                        style_function=lambda x: {'color':'transparent','fillColor':'transparent','weight':0},                        
                        tooltip=folium.features.GeoJsonTooltip(
                                      fields=['t.v','density'],
                                      aliases = ['區域村里','人口密度(人/公頃)'],
                                      labels=True,
                                      sticky=True  #顯示資訊方式不依樣(黏著度)
                                      )
                       ).add_to(m)


HeatMap(locations,name='HeatMap Covid-19 Case',radius=30).add_to(m)    #假設熱點半徑為30km


# 加入圖層開關功能 layer control 
folium.LayerControl().add_to(m)

m

Make this Notebook Trusted to load map: File -> Trust Notebook

outfp = "E:\Py_NTU 341\Homework\covid19heatmap&pop_density.html"
m.save(outfp)